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Abstract. We present an axisymmetric numerical model of a dynamo active accretion disc. If the dynamo-generated magnetic 
field in the disc is sufficiently strong (close to equipartition with thermal energy), a fast magneto-centrifugally driven outflow 
develops within a conical shell near the rotation axis, together with a slower pressure driven outflow from the outer parts of the 
disc as well as around the axis. Our results show that a dynamo active accretion disc can contribute to driving an outflow even 
without any external magnetic field. The fast outflow in the conical shell is confined by the azimuthal field which is produced 
by the dynamo in the disc and advected to the disc corona. This part of the outflow has high angular momentum and is cooler 
and less dense than its surroundings. The conical shell's half-opening angle is typically about 30° near the disc and decreases 
slightly with height. The slow outflow is hotter and denser. 
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1. Introduction 

The importance of magnetic fields for disc accretion is widely 
recognized (e.g., Mestel 1999), and the turbulent dynamo is be- 
lieved to be a major source of magnetic fields in accretion discs 
(Pudritz 1981a, 1981b; Stepinski & Levy 1988; Brandenburg et 
al. 1995; Hawley et al. 1996). A magnetic field anchored in the 
disc is further considered to be a major factor in launching and 
coUimating a wind in young stellar objects and active galactic 
nuclei (Blandford & Payne 1982; Pelletier & Pudritz 1992); see 
Konigl & Pudritz (2000) for a recent review of stellar outflows. 
Yet, most models of the formation and coUimation of jets rely 
on an externally imposed poloidal magnetic field and disregard 
any field produced in the disc. Our aim here is to study outflows 
in connection with dynamo-generated magnetic fields. We dis- 
cuss parameters of young stellar objects in our estimates, but 
the model also applies to systems containing a compact central 
object after rescaling and possibly other modifications such as 
the appropriate choice of the gravitational potential. 

Extensive numerical studies of collimated disc winds have 
been performed using several types of model. Uchida & 
Shibata (1985), Matsumoto et al. (1996) and Kudoh et al. 
(1998) consider an ideal MHD model of an accretion disc 
embedded in a non-rotating corona, permeated by an exter- 
nal magnetic field initially aligned with the disc rotation axis. 
Intense accretion develops in the disc (Stone & Norman 1994), 
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accompanied by a centrifugally driven outflow. The wind is 
eventually collimated by toroidal magnetic field produced in 
the corona by winding up the poloidal field. 

Bell (1994), Bell & Lucek (1995) and Lucek & BeU (1996, 
1997) use two- and three-dimensional ideal MHD models, with 
a polytropic equation of state, to study the formation and sta- 
bility of pressure driven jets collimated by poloidal magnetic 
field, which has a minimum inside the jet. The general struc- 
ture of thermally driven disc winds was studied analytically by, 
e.g., Fukue (1989) and Takahara et al. (1989). 

Another type of model was developed by Ustyugova et al. 
(1995), Romanova et al. (1997, 1998), Ouyed et al. (1997) 
and Ouyed & Pudritz (1997a, 1997b, 1999) who consider ideal 
MHD in a polytropic corona permeated by an external poloidal 
magnetic field and subsume the physics of the accretion disc 
into boundary conditions at the base of the corona. The disc 
is assumed to be in Keplerian rotation (with any accretion ne- 
glected). The corona is non-rotating initially, and the system 
is driven by the injection of matter through the boundary that 
represents the disc surface. These models develop a steady (or 
at least statistically steady) state consistent with the analyti- 
cal models by Blandford & Payne (1982), showing a magneto- 
centrifugal wind collimated by toroidal magnetic field, which 
again is produced in the corona by the vertical shear in the an- 
gular velocity. Three-dimensional simulations suggest that the 
resulting collimated outflow does not break due to the kink in- 
stability (Ouyed et al. 2003; Thomsen & Nordlund 2003). 

It is not quite clear how strong the external magnetic field 
in accretion discs can be. Dragging of an external field from 
large radii in a viscous disc requires that magnetic diffusiv- 
ity is much smaller than the kinetic one (Lubow et al. 1994), 
i.e. that the magnetic Prandtl number is significantly larger 
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than unity which would be difficuh to explain in a turbulent 
disc (Heyvaerts et al. 1996). This argument neglects, however, 
the effect of magnetic torques which could produce significant 
field line dragging even when the magnetic Prandtl number is 
of order unity (Shalybkov & Riidiger 2000). On the other hand, 
the efficiency of trapping an external magnetic field at initial 
stages of the disc formation is questionable because only a 
small fraction of the external flux can be retained if the den- 
sity contrast between the disc and the surrounding medium 
is large. In that case the magnetic field will be strongly bent 
and reconnection will remove most of the flux from the disc. 
Furthermore, the ionization fraction is probably too small in 
the inner disc to ensure sufficient coupling between the gas and 
flie magnetic field (Fromang et al. 2002). For these reasons it 
seems appropriate to explore whether or not magnetic fields 
generated by a dynamo within the disc can produce an outflow 
with realistic properties. 

The models of Ustyugova et al. (1995) and Ouyed & 
Pudritz (1997a, 1997b, 1999) have a distributed mass source 
and imposed poloidal velocity at the (unresolved) disc sur- 
face. As a result, their system develops a persistent outflow 
which is not just a transient. On the other hand, the models 
of Matsumoto et al. (1996) resolve the disc, but have no mass 
replenishment to compensate losses via the outflow and accre- 
tion, and so the disc disappears with time. Our model is an at- 
tempt to combine advantages of both of these approaches and 
also to add dynamo action in the disc. Instead of a rigidly pre- 
scribed mass injection, we allow for self-regulatory replenish- 
ment of matter within a resolved disc. Instead of prescribing 
poloidal velocity at the disc surface, we resolve the disc and 
prescribe an entropy contrast between the disc and the corona, 
leaving more freedom for the system. Such an increase of en- 
tropy with height is only natural to expect for a disc surrounded 
by a hot corona, and we parameterize the coronal heating by a 
(fixed) entropy contrast. We further add self-sustained, intrin- 
sic magnetic field to our system, as opposed to an external field 
used in the other models. Since our model goes beyond ideal 
MHD, the magnetic field must be maintained against decay. A 
simple form of mean-field dynamo action is included for this 
purpose. 

Like in the model of Ouyed & Pudritz (1997a, 1997b, 
1999), the hot, pressure supported corona does not rotate ini- 
tially. The disc is cool and is therefore centrifugally sup- 
ported, so its rotation is nearly Keplerian. The corresponding 
steady-state solution is used as our initial condition. This so- 
lution is, however, unstable because of the vertical shear in 
the angular velocity between the disc and the corona (Urpin 
& Brandenburg 1998). Angular momentum transfer by viscous 
and magnetic stresses also leads to a departure from the initial 
state. As a result, a meridional flow develops, which exchanges 
matter between the corona and the disc surface layers. Mass 
losses through the disc surface and to the accreting central ob- 
ject are then replenished in the disc where we allow local mass 
production whenever and wherever the density decreases be- 
low a reference value. Matter is heated as it moves into the 
corona; this leads to an increase in pressure which drives the 
wind. Another efficient driver of the outflow in our model is 
magneto-centrifugal acceleration. A strong toroidal magnetic 
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Fig. 1. General structure of the outflows typically obtained in 
the magnetic runs with mass sink presented in this paper. The 
cool, dense disc emits (a) a hardly collimated, thermally driven 
wind (slow, hot, dense, magnetized, rotating), (b) a magneto- 
centrifugally driven outflow in a conical shell (faster, cooler, 
less dense, magnetized, quickly rotating), and (c) a thermally 
driven outflow near the axis (slow, hot, dense, weakly magne- 
tized and weakly rotating). 



field produced by the dynamo in the disc is advected into the 
corona and contributes to confining the wind. 

Altogether, our model contains many features included in 
a range of earlier models. For example, the disc is essentially 
resistive to admit dynamo action, whereas magnetic diffusion 
turns out to be relatively unimportant in the corona, similarly 
to the models of, e.g., Wardle & Konigl (1993), Ferreira & 
Pefletier (1995) and Casse & FeiTeira (2000a, 2000b). The syn- 
thetic nature of the outflow, driven by both pressure forces and 
magneto-centrifugally, is characteristic of the self-similar so- 
lutions of Ferreira & Pelletier (1995) and Casse & Ferreira 
(2000b). The latter authors also stress the role of the verti- 
cal entropy gradient in enhancing the mass flux in the outflow 
and prescribe a vertical profile of the entropy production rate. 
Structured outflows outwardly similar to our results have been 
discussed by Krasnopolsky et al. (1999) and Goodson et al. 
(1997). However, in the former paper the structure in their so- 
lutions is due to the boundary conditions imposed on the disc 
surface, and in the latter model the structure is due to the inter- 
action between the stellar magnetic field and the disc under the 
presence of differential rotation between star and disc. On the 
contrary, the structure in our model results from the dominance 
of different driving mechanisms in the inner and outer parts of 
the disc. 

While we had initially expected to find collimated outflows 
similar to those reported by Ouyed & Pudritz (1997a, 1997b), 
the outflow patterns we obtained turned out to be of quite dif- 
ferent nature. However, a significant difference in our model is 
the resolved disc, which is necessary to model dynamo action 
in the disc, and the relatively small extension into the corona, 
where collimation is not yet expected. Figure [l| shows a sketch 
of the overall, multi-component structure of the outflows in the 
presence of magnetic fields and mass sink, which needs to be 
compared to the individual Figures in the rest of the paper. 
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The plan of the paper is as follows. We introduce our model 
in Sect. ^ and then consider a range of parameters in Sect. |] 
to clarify and illustrate the main physical effects. In Sect. ^ we 
explore the parameter space, and our conclusions are presented 
in Sect. || 

2. The model 

2. 1. Basic equations 

Our intention is to make our model as simple as possible and 
to avoid detailed modelling of mass supply to the accretion 
disc, which occurs differently in different accreting systems. 
For example, matter enters the accretion disc at large radii in 
a restricted range of azimuthal angles in binary systems with 
Roche lobe overflow. On the other hand, matter supply can be 
more uniform in both azimuth and radius in active galactic nu- 
clei. Being interested in other aspects of accretion physics, we 
prefer to avoid detailed modelling of these processes. Instead, 
similarly to the model of Ouyed & Pudritz (1997a, 1997b), we 
inject matter into the system, but an important difference is that 
we introduce a self-regulating mass source in the disc. We also 
allow for a mass sink at the centre to model accretion onto the 
central object (star). With these two effects included, the conti- 
nuity equation becomes 



9q „ 



(gu) 



(1) 



where u is the velocity field and g is the gas density. The mass 
source, Qg^^"^, is localized in the disc (excluding the star) and 
turns on once the local density drops below a reference value 
go (chosen to be equal to the initial density, see Sect. 2.4), i.e. 



disc ^disc('') Cstar('") / n 

Qg = (£'0 - + , 

'^disc 



(2) 



where the subscript '+' indicates that only positive values are 
used, i.e. 



x+ 



X if a; > 0, 
otherwise. 



(3) 



This means that matter is injected only if g < go, and that the 
strength of the mass source is proportional to the gas density 
deficit. Throughout this work we use cylindrical polar coordi- 
nates (tn, (p, z), assuming axisymmetry. The shapes of the disc 
and the central object are specified with the profiles ^disc and 

Cstar via 



Cdisc(K7, z) = 9(tuo-tn, d) Q{zo~\z\,d). 



(4) 
(5) 



Here, Q{x,d) is a smoothed Heaviside step function with a 
smoothing half- width d set to 8 grid zones in f disc and d — tq 
in Cstar; r — + is the spherical radius and ro is the 

stellar radius introduced in Eq. (^3]); -cuq and zq are the disc 
outer radius and disc half-thickness, respectively. So, ^disc is 
equal to unity at the disc midplane and vanishes in the corona, 
whereas ^stai is unity at the centre of the central object and van- 
ishes outside the star Note that ^disc — ^star > everywhere, 
since ^disc > 0) all the way to the origin, z = w = 0. 



In Eq. (Q), Tdisc is a response time which is chosen to be sig- 
nificantly shorter than the time scale of the depletion processes, 
which is equivalent to the time scale of mass replenishment 
in the disc. A/disc /Msourcc (cf. Sect. 3.4), to avoid unphysical 
influences of the mass source. We do not fix the distribution 
and magnitude of g^'^'^ beforehand, but the system adjusts it- 
self such as to prevent the disc from disappearing. 

The self-regulating mass sink at the position of the central 
star is modelled in a similar manner. 



Cstar(T-) 



{q- Qo)- 



(6) 



where ^star is defined in Eq. (g) and Tgtar is a central accretion 
time scale that controls the efficiency of the sink. We discuss 
physically meaningful values of Tgtai in Sect. 2.5. 

Apart from the continuity equation (|l]), the mass source 
also appears in the Navier-Stokes equation, unless matter is 
always injected with the ambient velocity of the gas. In that 
case, however, a runaway instability can occur: if matter is al- 
ready slower than Keplerian, it falls inward, and so does the 
newly injected matter. This enhances the need for even more 
mass injection. A similar argument applies also if matter is ro- 
tating faster than Keplerian. This is why we inject matter at 
the Keplerian speed. This leads to an extra term in the Navier- 
Stokes equation, (mk — m)(7^"''^, which would only be absent if 
the gas were rotating at the Keplerian speed. Thus, the Navier- 
Stokes equation takes the form 

^ = -ivp-V$ + i[F + («K-«)gf1, (7) 

where D/Dt — d/dt+u-'V is the advective derivative, p is the 
gas pressure, $ is the gravitational potential, = JxB+'V t 
is the sum of the Lorentz and viscous forces, J = V x i?//io 
is the current density due to the (mean) magnetic field B, /io is 
the magnetic permeability, and x is the viscous stress tensor. 



gv{duil dxj + duj/dxi). 



(8) 



Here, v is the kinematic viscosity, which has been subdivided 
into three contributions. 



1^ = 1^1+ i^adv + i^shock- 



(9) 



The first term is a turbulent (Shakura-Sunyaev) viscosity in the 
disc. 



i^t = assCsZo^discir), 



(10) 



where Cg = (773/ gY^^ is the sound speed and 7 is the ratio of 
specific heats. The second term is an artificial advection vis- 
cosity. 



(11) 



which is required to stabilize rapidly moving patterns. Here, 

6x = min((5-n7, 6z) is the mesh size, and c^*^^ is a con- 



stant specified in Sect. 2.5 In Eq. (njl) we have only used 



the poloidal velocity, Mpoi, and the Alfven speed due to the 
poloidal magnetic field, UA,poi — [Bp^i/ g^oY^^, because in 
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an axisymmetric calculation advection in the cp-direction is 
unimportant. Finally, 



^shock 



^shock 



5x' 



(12) 



is an artificial shock viscosity, with a constant specified 

in Sect. 2.5. Note that v^^v and i^shock are needed for numerical 
reasons; they tend to zero for increasing resolution, 5x 0. 

We assume that the magnetic field in the disc is gener- 
ated by a standard a^fi-dynamo (e.g., Krause & Radler 1980), 
which implies an extra electromotive force, aB, in the induc- 
tion equation for the mean magnetic field, B. To ensure that B 
is solenoidal, we solve the induction equation in terms of the 
vector potential A, 



dA 
'dt 



= u X B + aB — ti^qJ, 



(13) 



where B — x A, and 77 is the magnetic diffusivity. 

Since a{r) has to be antisymmetric about the midplane and 
vanishing outside the disc, we adopt the form 



ao 



Z_ Cdisc(^) 

zq 1 + vl/v^ 



(14) 



where wa is the Alfven speed based on the total magnetic field, 
and q;o and vq are parameters that control the intensity of dy- 
namo action and the field strength in the disc, respectively. The 
a-effect has been truncated near the axis, so that a = for 
w < 0.2. For the magnetic diffusivity we assume 



V = Vo +Vt, 



(15) 



where rjo is a uniform background diffusivity and the turbulent 
part, rjt = '7to$disc('"), vanishes outside the disc. Thus, mag- 
netic diffusivity in the corona is smaller than in the disc. 

Depending on the sign of a and the vertical distribution of 
77, the dynamo can generate magnetic fields of either dipolar or 
quadrupolar symmetry. We shall discuss both types of geome- 
try. 

2.2. Implementation of a cool disc embedded in a 
hot corona 

Protostellar systems are known to be strong X-ray sources 
(see, e.g., Glassgold et al. 2000; Feigelson & Montmerle 1999; 
Grosso et al. 2000). The X-ray emission is generally attributed 
to coronae of disc-star systems, plausibly heated by small scale 
magnetic reconnection events (Galeev et al. 1979), for exam- 
ple in the form of nanoflares that are caused by slow footpoint 
motions (Parker 1994). Heating of disc coronae by fluctuating 
magnetic fields is indeed quite natural if one accepts that the 
disc turbulence is caused by the magneto-rotational instabil- 
ity. Estimates for the coronal temperatures of YSOs range from 
IQS K to 10^ K (see, e.g., Feigelson & Montmerle 1999). For 
the base of the disc corona, temperatures of at least 8 x 10"^ 
are to be expected in order to explain the observed mass loss 
rates (Kwan & Tademaru 1995; Kwan 1997). The discs, on 
the other hand, have typical temperatures of a few 10"^ K (e.g., 
Papaloizou & Terquem 1999). 



A simple way to implement a dense, relatively cool disc 
embedded in a rarefied, hot corona without modelling the de- 
tailed physics of coronal heating is to prescribe the distribution 
of specific entropy, s(r), such that s is smaller within the disc 
and larger in the corona. For a perfect gas this implies p = Kg'* 
(in a dimensionless form), where K = e^l'^^ is a function of 
position (here p and q are gas pressure and density, 7 = Cp/cy, 
and Cp and are the specific heats at constant pressure and 
constant volume, respectively). 

We prescribe the polytrope parameter K to be unity in the 
corona and smaller in the disc, so we put 

[K{r)Yh = e^/cp ^ 1 _ (l-/3)edi,c(r), (16) 

where < < 1 is a free parameter with (3 — e^'^^/^p, 
that controls the entropy contrast. As > 0, between corona 
and disc. We consider values of (3 between 0.1 and 0.005, 
which yields an entropy contrast, As/cp, between 2.3 and 5.3. 
The temperature ratio between disc and corona is roughly (3; 
see Eq. (26|). Assuming pressure equilibrium between disc and 
corona, and p (x pT for a perfect gas, the corresponding den- 
sity ratio is f3^^. 

2.3. Formulation in terms of potential enthalpy 

In the present case it is advantageous to use the potential en- 
thalpy. 



(17) 



as a variable. Here, h is the specific enthalpy, h = c^T+p/ g = 
CpT for a perfect gas (with constant specific heats), and T is 
temperature. Therefore, specific enthalpy h is related to p and 
g via h = 7(7— l)^^p/ g- Specific entropy s is related to p and 
g (up to an additive constant) through s — Cy\w-p ~ Cp\n g for 
a perfect gas. We have therefore 

Dln/i Dlnp D In (? , ^D\ng 7 Ds 

(7-I)— ft;— + — — ■ (18) 



Dt 



Dt 



Dt 



Dt 



Cp Dt 



Since sis independent of time, Ds/Dt = it- Vs. Together with 
Eq. (|l]), this yields an evolution equation for h, which can be 
written in terms of H as 



= u ■ V$ + (7 - ( ^ - V • M 
Dt \g 



7/1 



u-Vs. (19) 



In the following we solve Eq. (|T9|) instead of Eq. (|l]). In terms 
of h, density and sound speed are given by 

.7-1 ^ Tzl A , 

7 K' 



Q 



(7 - = 7 



P 



(20) 



It proved necessary to include an artificial diffusion term in 
Eq. ( [T^ ) with a diffusion coefficient proportional to u. 

The first law of thermodynamics allows us to express the 
pressure gradient in terms of h and s, 

--Vp=-Vh + TVs, (21) 
Q 

and with T = h/cpWe obtain 

^ = - Vi7 + hVs/cp + i [F + («K - u)qf'-] , (22) 
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which now replaces Eq. (^. 

We use a softened, spherically symmetric gravitational po- 
tential of the form 



-GM, (r" + r^) 



,n\-l/n 



(23) 



where G is the gravitational constant, is the mass of the 
central object, r is the spherical radius, tq is the softening ra- 
dius, and n = 5; tentatively, ro can be identified with the stellar 
radius. 



2.4. The initial state 

Our initial state is the hydrostatic equilibrium obtained by solv- 
ing, for h, the vertical balance equation obtained from Eq. (0), 



oz 



oz 



(24) 



from large z (where /i + $ = 0) down to z = 0. The initial den- 
sity distribution go is then obtained using Eq. (pO[); it decreases 
monotonically with both z and w in the equilibrium state. 

The initial rotation velocity, u^p^, follows from the radial 
balance equation. 



ds/cp 
dm 



(25) 



In the disc, h — CpT is small, so u^o is close to the Keplerian 
velocity, while the corona does not rotate initially, and so is 
supported by the pressure gradient. 

As a rough estimate, the value of h in the midplane of the 
disc is 

/idisc ~ ^ ^corona = /3GM*/a7, (26) 

as can be seen by integrating Eq. (Q), ignoring the d^/dz 
term. Thus, the initial toroidal velocity in the midplane can 
be obtained from Eq. ( |2^ ) using Eq. and recalling that 
ds /dm = in the midplane, which gives 



Vo ~ v/(l - l3)GM,/m = VW?WK. 



(27) 



For /? = 0.1, for example, the toroidal velocity is within 5% of 
the Keplerian velocity. 

2.5. Dimensionless variables and choice of 
parameters 

Our model is scale invariant and can therefore be applied to 
various astrophysical objects. We consider here the range of 
parameters typical of protostellar discs, for which a typical sur- 
face density is Sq ~ 1 gcm~^. A typical coronal sound speed 
is Cso ~ lO^kms^^, which corresponds to a temperature of 
r « 4 X 10^ K. This allows us to fix relevant units as follows: 



[S] = So, M = Cso, [x] = GAU 



(28) 



With M, = IMq we have [x] « 0.1 AU, and [t] = 
« 1.5 day. The unit for density is [g] — \n\/[x\ « 
7 X lO^^'^ g cm^"^, the unit for the mass accretion rate is 
[M] = « 2 x lO-'^Mgyr^i, and the unit for the 

magnetic field is [B] = [u]{fio[g])'^^'^ « 30 G. 



Since [h] = [u]^, the dimensionless value h — 1 coiTe- 
sponds to 10^^ cm^ s^^. With a mean specific weight fi — 0.6, 
the universal gas constant TZ = 8.3 x lO^cm^s^^K^^ and 
7 = 5/3, we have 



7 7^ 

7- 1 ^ 



3.5 X lO^cm^s-^K-^ 



(29) 



Therefore, h — 1 corresponds to T = [u]^/[s] = [h]/cp w 
3 X 10^ K. Using h = in the corona, this corresponds to a 
temperature of 3 x 10^ K at r = [x] = 0.1 AU. 

We choose /3 between 0.1 and 0.005, corresponding to 
a typical disc temperature (in the model) of 3 x 10"* K to 
1.5 X 10'^ K; see Eq. (p6|); real protostellar discs have typical 
temperatures of about a few thousand Kelvin (see Sect. 

In our models, we use mo = 1.5 for the disc outer radius, 
zo = 0.15 for the disc semi-thickness and ro = 0.05 for the 
softening (stellar) radius. The disc aspect ratio is zo/mo = 0.1. 
Note that ro = 0.05 coiTesponds to 7 x 10^*^ cm, i.e. one solar 
radius. Therefore, we shall not reduce it much below this phys- 
ically meaningful value. Note, however, that smaller values of 
ro would result in faster outflows (Ouyed & Pudritz 1999). 

Furthermore, cf"" = 0.02 and 4'^°'='^ = 1.2. We vary the 
value of ass between 0.003 and 0.007. 

The mean-field dynamo is characterized by the parameters 
laol = 0.3, vq — Cs, rjta — 10^^ and t^o = 5 x lO^"'. The 
total magnetic diffusivity in the disc is therefore rjTo = jyto + 
rjQ = 0.0015. In terms of the usual Shakura-Sunyaev viscosity 
parameter, this coiTesponds to 



(77) _ VTO 



^SS 



0.01 f^V^' f ^^u^r ; (30) 

V 13 J V0.0015/ V0.15/ ' ^ ' 



""s.disc 

-1/2' 



s, corona 



[cf. Eqs @ and @] 



Cs,disc^O 

where we have used c? 

and Cg^ corona ~ ^ 

In terms of the usual nondimensional dynamo parameters 
we have 



\Ca \ = \ao\zo/r]To = 30 

and, for Keplerian rotation, 

G^ = imdn/dm)z^/rjTo = -22.5m-^/'^, 

so that the dynamo number is given by 

\V\ = = 675tn-3/2. 



(31) 



(32) 



(33) 



Note that the value of the dynamo number expected for accre- 
tion discs is given by 



\V\ = 



\ao mdft/dm\ Zq 



4m 



(34) 



for ?7to = alJ'^CsAiscZo,\ao\ c^^disc and Csaisc = ^zq. As 
follows from Eqs (|o|)-( 33 ), for Keplerian rotation. 



\V 



[4iy ^\ao\zop-'m-'^'. 



(35) 



[Note that this expression is independent of 77x0-] For our 
choice of parameters of |q;o| — 0.3 and zq = 0.15, expres- 
sion @ is < 1 if ti7 ^ 0.002 /3-2 = 0.2 for /? = 0.1, which 
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corresponds to the truncation radius of the a-effect. Therefore, 
our choice of parameters is consistent with the constraint (p4|). 

We take 7 = 5/3, Tdisc = 0.1, and consider two values of 
Tstar- For Tstar oo, the mass sink at the central object is sup- 
pressed, whereas Tstar implies instantaneous accretion of 
any extra matter (relative to the hydrostatic equilibrium) by the 
central object. A realistic lower limit for Tstar can be estimated 
as Tstar ^ To/vff w 0.008, whcrc Wff is the free fall velocity 
(given by iw| = GM^^/tq in dimensional quantities). In most 
cases we used Tgtar = 0.01, but we also tried an even smaller 
value of Tstar = 0.005 and obtained very similar results. A fi- 
nite value of Tstar implies that matter is not instantaneously ab- 
sorbed by the sink. Therefore, some matter can leave the sink 
if it moves so fast that its residence time in the sink is shorter 
than Tstar. As can be seen below, a small (negligible) fraction 
of mass does indeed escape from the sink. 

Computations have been carried out in domains ranging 
from {m, z) e [0, 2] x [-1, 1] to [0, 8] x [-2, 30] , but the results 
are hardly different in the overlapping parts of the domains. In 
our standard computational box, 5w = 5z = 0.01 and in the 
case of the larger computational domain 5w = 5z = 0.02. 




-1„.9 -Q.5 a.Q Q.5 1.0 



2.6. Numerical method and boundary conditions 

We use third order Runge-Kutta time-stepping and a sixth or- 
der finite-difference scheme in space. Details and test calcula- 
tions are discussed by Brandenburg (2001). 

On the outer boundaries, the induction equation is evolved 
using one-sided derivatives (open boundary conditions). The 
normal velocity component has zero derivative normal to the 
boundary, but the velocity is required to be always directed 
outwards. The tangential velocity components and potential en- 
thalpy on the boundaries are similarly obtained from the next 
two interior points away from the boundary. Tests show that 
the presence of the boundaries does not affect the flow inside 
the computational domain. Regularity conditions are adopted 
on the axis where the zu and components of vectorial quan- 
tities vanish, whereas scalar variables and the z components of 
vectorial quantities have vanishing radial derivative. 



3. Results 

In this section we discuss a range of models of increasing de- 
gree of complexity. We first consider in Sect. 3.1 the simplest 
model that contains neither magnetic field nor mass sink at the 
centre to show how the pressure gradient resulting from the 
(fixed) entropy distribution drives a disc outflow. It is further 
shown that the outflow is significantly reduced if accretion onto 
the central object is allowed, but restored again if the disc can 
generate a large-scale magnetic field (Sect. 3.2). Having thus 
demonstrated the importance of the large-scale magnetic field 
in our model, we discuss in Sect. 3.3 its structure, largely con- 
trolled by the dynamo action but affected by the outflow. Model 
parameters used in these sections are not necessarily realistic as 
we aim to illustrate the general physical nature of our solutions. 
We present a physically motivated model in Sect. ^ where the 
set of model parameters is close to that of a standard accretion 



Fig. 2. A nonmagnetic model without mass sink at the cen- 
tre: velocity vectors and logarithmically spaced density con- 
tours (from 0.008 to 6000, separation factor of 1.6) super- 
imposed on a grey scale representation of logip h at time 
t = 223. Specific enthalpy h is directly proportional to tem- 
perature T, and logj^g /i = (—1,0,1) corresponds to T « 
(3x10'*, 3x10^ 3x10*^) K. Note hot gas near the central ob- 
ject and in the near corona, and cooler gas in the disc and the 
far corona. The dashed line shows the sonic surface where the 
poloidal velocity equals Cs. The disc boundary is shown with a 
thin black line; /3 = 0.1. 



disc around a protostellar object. The physical nature of our 



solutions is discussed in Sects. 3.4, 3.5 and 



3. 1. Nonmagnetic outflows 

We illustrate in Fig. ^ results obtained for a model without any 
magnetic fields and without mass sink. A strong outflow de- 
velops even in this case, which is driven mostly by the vertical 
pressure gradient in the transition layer between the disc and 
the corona, in particular by the term TVs in Eq. (|T]). The gain 
in velocity is controlled by the total specific entropy difference 
between the disc and the corona, but not by the thickness d of 
the transition layer in the disc profile (Q). 

The flow is fastest along the rotation axis and within a cone 
of polar angle of about 30°, where the terminal velocity « 3 
is reached. The conical shape of the outflow is partly due to 
obstruction from the dense disc, making it easier for matter to 
leave the disc near the axis. Both temperature and density in 
nonmagnetic runs without mass sink are reduced very close to 
the axis where the flow speed is highest. 

The general flow pattern is sensitive to whether or not mat- 
ter can accrete onto the central object. We show in Fig. || re- 
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Fig. 3. As in Fig. ||, but with a mass sink at the centre, 
Tatar = 0.01, at time t — 312. The outflow speed near the axis 
is strongly reduced in comparison to that in the model without 
mass sink shown in Fig. ^. 

suits for the same model as in Fig. ^ but with a mass sink 
given by Eq. (^) with Tgtar = 0.01. This can be compared with 
earlier work on thermally driven winds by Fukue (1989) who 
also considered polytropic outflows, but the disc was treated 
as a boundary condition. In Fukue (1989), outflows are driven 
when the injected energy is above a critical value. The origin of 
this energy injection may be a hot corona. The critical surface 
in his model (see the lower panel of his Fig. 2) is quite similar 
to that found in our simulation (our Fig. |]), although our open- 
ing angle was found to be larger than in Fukue's (1989) model. 
[Below we show, albeit with magnetic fields, that smaller val- 
ues of P do result in smaller opening angles, see Fig. ^ which 
would then be compatible with the result of Fukue (1989).] 

As could be expected, the mass sink hampers the outflow 
in the cone (but not at tu ^ 0.5). The flow remains very sim- 
ilar to that of Fig. ^ when Tstar is reduced to 0.005. Thus, the 
nonmagnetic outflows are very sensitive to the presence of the 
central sink. As we show now, magnetized outflows are differ- 
ent in this respect. 

3.2. Magnetized outflows 

In this section we discuss results obtained with magnetic fields, 
first without mass sink at the centre and then including a sink. 
We show that the effects of the sink are significantly weaker 
than in the nonmagnetic case. 

A magnetized outflow without the central mass sink, shown 
in Fig. Q is similar to that in Fig. |^, but is denser and hotter 
near the axis, and the high speed cone has a somewhat larger 
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Fig. 4. A magnetic model without mass sink at the centre: ve- 
locity vectors and poloidal magnetic field lines (white) super- 
imposed on a grey scale representation of temperature (in terms 
of h) for a run with ao = —0.3 (resulting in roughly dipolar 
magnetic symmetry) at time t — 269. In contrast to the non- 
magnetic model of Fig. ^ a conical shell has developed that is 
cooler and less dense than its exterior. The conical shell inter- 
sects z = ±1 at around w « 1.2. The dashed line shows the 
surface where the poloidal velocity equals (cg^ + v\ poi)^^^^ 
with t^A.poi the Alfven speed from the poloidal magnetic field 
(fast magnetosonic surface with respect to the poloidal field). 
The disc boundary is shown with a thin black line; (3 — 0.1. 



opening angle. In addition, the outflow becomes more struc- 
tured, with a well pronounced conical shell where tempera- 
ture and density are smaller than elsewhere (the conical shell 
reaches z = ±1 at sa 1.2 in Fig. Here and in some of the 
following figures we also show the fast magnetosonic surface 
with respect to the poloidal field. In Sect. we show that this 
surface is close to the fast magnetosonic surface. As shown in 
Fig. ^ the outflow becomes faster inside the cone (uz ~ 5 on 
the axis). As expected, we find that deeper potential wells, i.e. 
smaller values of tq in Eqs (jsj) and (^), result in even faster 
flows and in larger opening angles. 

Our results are insensitive to the size and symmetry of the 
computational domain: we illustrate this in Fig. ^ with a larger 
domain of the size [0, 8] x [—2, 30]. The disc midplane in this 
run is located asymmetrically in z in order to verify that the 
(approximate) symmetry of the solutions is not imposed by the 
symmetry of the computational domain. Figure || confirms that 
our results are not affected by what happens near the computa- 
tional boundaries. 

Unlike the nonmagnetized system, the magnetized outflow 
changes only comparatively little when the mass sink is intro- 
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Fig. 6. Results for a larger domain, [0, 8] x [—2, 30], with the same parameters as in Fig. ^ averaged over times t = 200 . . . 230, 
P = 0.1. Left panel: velocity vectors, poloidal magnetic field lines and grey scale representation of h in the inner part of the 
domain. Right panel; velocity vectors, poloidal magnetic field lines and normalized specific enthalpy h/\^\ in the full domain. 



duced at the centre. We show in Fig. ^ the results with a sink 
(TVitar — 0.01) and otherwise the same parameters as in Fig. ^. 
As could be expected, the sink leads to a reduction in the out- 
flow speed near the axis; the flow in the high speed cone be- 
comes slower But apart from that, the most important effects 
of the sink are the enhancement of the conical structure of the 
outflow and the smaller opening angle of the conical shell. A 
decrease in Tstai by a factor of 2 to 0.005 has very little effect, 
as illustrated in Figs. p|andp^ 



Increasing the entropy contrast (while keeping the specific 
entropy unchanged in the corona) reduces the opening angle of 
the conical shell. Pressure driving is obviously more important 
in this case, as compared to magneto-centrifugal driving (see 
Sect. 3.5). A model with /3 — 0.02 (corresponding to a density 
and inverse temperature contrast of about 50: 1 between the disc 
and the corona) is shown in Fig. ^ At i9 = 60°, the radial 
velocity Ur is slightly enhanced relative to the case /3 = 0.1 
(contrast 10:1); see Fig. |l^ At = 30°, on the other hand, the 
flow with the larger entropy contrast reaches the Alfven point 



B. von Rekowski et al.; Structured outflow from a dynamo active accretion disc 



9 



r(,=0.02, magn., ao<0 
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Fig. 5. Profiles of vertical velocity and sound speed along the 
axis, = 0, for the nonmagnetic model of Fig. ^ with smooth- 
ing radius ro = 0.05 (dotted) and two runs with magnetic field, 
with ro = 0.05 as in Fig. ^ (solid) and ro = 0.02 (dashed, 
t = 162). The shaded area marks the location of the disc. 
Terminal wind speeds are reached after approximately three 
disc heights. The presence of a magnetic field and a deeper po- 
tential well (smaller rp) both result in faster outflows. The mod- 
els shown here have no mass sink at the centre and P = 0.1. 



close to the disc (at r w 0.27 as opposed to r ^ 1 in the other 
case), which leads to a smaller terminal velocity. 

We conclude that the general structure of the magnetized 
flow and its typical parameters remain largely unaffected by 
the sink, provided its efficiency t^^^. does not exceed a certain 
threshold. It is plausibly the build-up of magnetic pressure at 
the centre that shields the central object to make the central ac- 
cretion inefficient. This shielding would be even stronger if we 
included a magnetosphere of the central object. We discuss in 
Sect. p~ 



the dependence of our solution on the geometrical size 
of the sink and show that the general structure of the outflow 
persists as long as the size of the sink does not exceed the disk 
thickness. 



3.3. Magnetic field structure 

The dynamo in most of our models has uq < 0, consis- 
tent with results from simulations of disc turbulence driven 
by the magneto-rotational instability (Brandenburg et al. 1995; 
Ziegler & Riidiger 2000). The resulting field symmetry is 
roughly dipolar, which seems to be typical of ail disc dynamos 
with ao < in a conducting corona (e.g., Brandenburg et al. 
1990). We note that the dominant parity of the magnetic field is 
sensitive to the magnetic diffusivity in the corona: a quadrupo- 
lar oscillatory magnetic field dominates for ao < if the disc 
is surrounded by vacuum (Stepinski & Levy 1988). 

For ao < 0, the critical value of ao for dynamo action 
is about 0.2, which is a factor of about 50 larger than without 
outflows. Our dynamo is then only less than twice supercritical. 
A survey of the dynamo regimes for similar models is given by 
Bardou et al. (2001). 

The initial magnetic field (poloidal, mixed parity) is weak 
[pinag = -B^/ (2//o) ~ lO^'^], cf. Fig. [Ill for comparison with 
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Fig. 7. As in Fig. ^ but with a mass sink at the centre, 
Tstar = 0.01. The outflow speed near the axis and the opening 
angle of the conical shell (now reaching z = ±1 at w 0.5) 
are reduced in comparison to that in the model without mass 
sink shown in Fig. ^ but the general structure of the outflow is 
little affected. Averaged over times t = 130 . . . 140, /3 = 0.1. 
[Reference model.] 



the gas pressure], but the dynamo soon amplifies the field in the 
disc to Pmag.tor ~ 10, and then supplies it to the corona. As a 
result, the corona is filled with a predominantly azimuthal field 
with Pmag,tor/p ~ 100 at larger radii; see Fig. [ri|. We note, 
however, that the flow in the corona varies significantly in both 
space and time]]. Magnetic pressure due to the toroidal field B^, 
Pmag,tor, excceds gas pressure in the corona outside the inner 
cone and confines the outflow in the conical shell. The main 
mechanisms producing in the corona are advection by the 
wind and magnetic buoyancy (cf. Moss et al. 1999). Magnetic 
diffusion and stretching of the poloidal field by vertical shear 
play a relatively unimportant role. 

The field in the inner parts of the disc is dominated by the 
toroidal component; \B^p/Bz\ ~ 3 at tu ^ 0.5; this ratio is 
larger in the corona at all ru. However, as shown in Fig. [ill, this 
ratio is closer to unity at larger radii in the disc. 

As expected, ap > results in mostly quadrupolar fields 
(e.g., Ruzmaikin et al. 1988). As shown in Fig. [T^, the mag- 
netic field in the corona is now mainly restricted to a narrow 
conical shell that crosses z = ±1 at to w 0.6. Comparing 
this figure with the results obtained with dipolar magnetic fields 
(Fig. one sees that the quadrupolar field has a weaker effect 
on the outflow than the dipolar field; the conical shell is less 



' See movie at http://www.nordita.dk/~brandenb/ 
movies/ outflow 
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Fig. 8. As in Fig. 0, but with a stronger mass sink at the centre, 
Tstar = 0.005. The flow pattern is very similar to that of Fig. ^ 
where the time scale of the sink is twice as large. Averaged over 
times t = 122... 132, /3 = 0.1. 
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Fig. 9. As in Fig. 0, but with a larger entropy contrast, (3 — 
0.02. The opening angle of the conical shell is reduced in com- 
parison to that of Figs. ^ and || where the entropy contrast is 
smaller (the shell crosses z = ±1 at tu sa 0.15). Averaged over 
times t = 140.. .240, r^tar = 0.01. 
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Fig. 10. Spherical radial velocity component, Ur (dashed), and 
poloidal Alfven speed, VA,po\ (solid), as functions of spherical 
radius at polar angles i? = 30° (thick lines) and = 60° (thin 
lines) for the models of Figs. ^ (top panel), ^ (middle panel) 
and ^ (bottom panel). 



pronounced. However, the structures within the inner cone are 
qualitatively similar to each other. 

The magnitude and distribution of a in Eq. ( |l4} ) only 
weakly affect magnetic field properties as far as the dynamo is 
saturated. For a saturated dynamo, the field distribution in the 
dynamo region (0.2 < m < 1.5, \z\ < 0.15) roughly follows 
from the equipartition field, B ~ {gjiov'^Y^^ with va = Cg- 
In other words, nonlinear states of disc dynamos are almost in- 
sensitive to the detailed properties of a (e.g.. Beck et al. 1996; 
Ruzmaikin et al. 1988). 

A discussion of disc dynamos with outflows, motivated by 
the present model, can be found in Bardou et al. (2001). It is 
shown that the value of magnetic diffusivity in the corona does 
not affect the dynamo solutions strongly. Moreover, the outflow 
is fast enough to have the magnetic Reynolds number in the 
corona larger than unity, which implies that ideal integrals of 



motion are very nearly constant along field lines; see Sect. 3.6 
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Fig. 11. Radial profiles of magnetic pressure from the toroidal 
(solid) and poloidal (dotted) magnetic fields and thermal pres- 
sure (dashed) for the model of Fig. ^. Shown are the averages 
over the disc volume (lower panel) and over a region of the 
same size around z = 8 in the corona (upper panel). /3 = 0.1. 
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Fig. 12. As in Fig. but with ao = +0.3, at time t = 132. The 
magnetic field geometry is now mostly quadrupolar because 
ao > 0. 



The most important property is the sign of a as it controls the 
global symmetry of the magnetic field. 



3.4. Mass and energy loss 

The mass injection and loss rates due to the source, sink and 
wind are defined as 



M^ourcc = / qf'' dV, Afsink = / 



and 

Mwind ^ (t QU- dS, 



(36) 



(37) 



respectively, where the integrals are taken over the full compu- 
tational domain or its boundary. About 1 /3 of the mass released 
goes into the wind and the rest is accreted by the sink, in the 
model with Tgtar = 0.01 and /3 — 0.1 of Fig. 0. Reducing Tgtar 
by a factor of 2 (as in the model of Fig. only changes the 
global accretion parameters by a negligible amount (^ 10%). 

The mass loss rate in the wind fluctuates on a time scale of 
5 time units, but remains constant on average at about A/wind ~ 
3, corresponding to 6 x 10~^ Mq yr~^, in the models of Figs. ^ 
and ^ The mass in the disc, Mdisc, remains roughly constant. 

The rate at which mass needs to be replenished in the disc, 
A/source A^disc, IS about 0.4. This rate is not controlled by the 
imposed response rate of the mass source, t^^^., which is 25 
times larger So, the mass source adjusts itself to the disc evo- 
lution and does not directly control the outflow. We show in 
Fig. |l3| trajectories that start in and around the mass injection 
region. The spatial distribution of the mass replenishment rate 
shown in Fig. [ij indicates that the mass is mainly injected 



^disc 



close to the mass sink, and q'^^^'^ remains moderate in the outer 
parts of the disc. (Note that the reduced effect of the mass sink 
in the magnetized flow is due to magnetic shielding rather than 
to mass replenishment near the sink - see Sect. 3.2.) 



The angular structure of the outflow can be characterized by 
the following quantities calculated for a particular spherical ra- 
dius, r — 8, for the model of Fig. ^: the azimuthally integrated 
normalized radial mass flux density, Af (?9)/A/disc, where 



M{-d) = 2TTr^gUrSm-d, Afdisc = / pdV, 

' disc 



(38) 



the azimuthally integrated normalized radial angular momen- 
tum flux density, J(i9) / Jdisc, where 

J[-d) = 2TTr'^ gvju^UrSind, Jdisc = / p-iuu^ dV, (39) 

J disc 

the azimuthally integrated normalized radial magnetic energy 
flux (Poynting flux) density, Eu_{'&) / Eu_, where 



■ , , -yiE X B)r 



Mo 



sin,?, Em^ / — dF, (40) 

disc ^MO 



and the azimuthally integrated normalized radial kinetic energy 
flux density, -Ek(i^) /-^k, where 

E]^{§) ^2TTr^{\QV?Ur)sm§, Ek ^ [ \pu^ dV. (41) 



Here, Afdisc, •/disc, and i?K are the mass, angular mo- 
mentum, and magnetic and kinetic energies in the disc. A polar 
diagram showing these distributions is presented in Fig. O for 
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Fig. 13. Azimuthally integrated mass flux density, represented as a vector 2Tr wq{utxj,Uz), in the simulation of Fig. |^ with a 
dipolar magnetic field and a mass sink at the centre. Shades of grey show the distribution of the mass source in the disc, Qg^'^'^. 
The disc boundary is shown with a white line. 



the model of Fig. ^. Note that this is a model without central 
mass sink where the flow in the hot, dense cone around the 
axis is fast. The fast flow in the hot, dense cone carries most of 
the mass and kinetic energy. A significant part of angular mo- 
mentum is carried away in the disc plane whilst the magnetic 
field is ejected at intermediate angles, where the conical shell 
is located. 

The radial kinetic and magnetic energy flux densities, in- 
tegrated over the whole sphere, are Lk = Jq Ek di? w 
0.54Mwindc2 and Lm = EuWdd w 0.03 AfwindC^, re- 
spectively, where « 2.8 is the fast magnetosonic speed (with 
respect to the poloidal field) at the critical surface (where Uz — 
c,) on the axis. Thus, M„indcl can be taken as a good indica- 
tor of the kinetic energy loss, and the magnetic energy loss into 
the exterior is about 3% of this value. These surface-integrated 
flux densities (or luminosities) are, as expected, roughly inde- 
pendent of the distance from the central object. 

3.5. Mechanisms of wind acceleration 

The magnetized outflows in our models with central mass sink 
have a well-pronounced structure, with a fast, cool and low- 
density flow in a conical shell, and a slower, hotter and denser 
flow near the axis and in the outer parts of the domain. Without 
central mass sink, there is a high speed, hot and dense cone 
around the axis. 

The magnetic field geometry (e.g.. Fig. 0) is such that for 
m > 0.1, the angle between the disc surface and the field 
lines is less than 60°, reaching w 20° at « 1-1.5, which is 
favourable for magneto-centrifugal acceleration (Blandford & 
Payne 1982; Campbell 1999, 2000, 2001). However, the Alfven 
surface is so close to the disc surface in the outer parts of the 
disc that acceleration there is mainly due to pressure gradient. 
The situation is, however, different in the conical shell contain- 
ing the fast wind. As can be seen from Fig. [T^, the Alfven sur- 
face is far away from the disc in that region and, on a given field 
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Fig. 14. Dependence, on polar angle of azimuthally inte- 
grated radial mass flux density M{d) through a sphere r — 
8 (solid, normalized by the disc mass Mdisc ~ 12), az- 
imuthally integrated radial angular momentum flux density 
J(??) (dashed-dotted, normalized by the disc angular momen- 
tum Jdisc ~ 4.3), azimuthally integrated radial Poynting flux 
density Emii^) (divided by 5, dashed, normalized by the mag- 
netic energy in the disc Em ~ 0.6), and azimuthally inte- 
grated radial kinetic energy flux density Ek^'^) (divided by 
10, dash-3dots, normalized by the kinetic energy in the disc 
Ek ~ 9.7), for the model of Fig. |[ The unit of all the quanti- 
ties is [t]~^ rad^^. 



line, the Alfven radius is at least a few times larger than the ra- 
dius of the footpoint in the disc. This is also seen in simulations 
of the magneto-centrifugally driven jets of Krasnopolsky et al. 
(1999); see their Fig. 4. The lever arm of about 3 is sufficient for 
magneto-centrifugal driving to dominate. As can be seen also 
from Fig. IT^ the flow at the polar angle -d = 60° is mainly ac- 
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Fig. 15. Angular momentum (normalized by the maximum 
angular momentum in the disc) for the model with Tstai ~ 
0.01 and /? = 0.1, shown in Fig. 0. The maximum value is 
V^discmax ~ 2.5. The dashed line shows the fast magnetosonic 
surface with respect to the poloidal field (cf. Fig. the solid 
line the Alfven surface where the poloidal velocity equals the 
poloidal Alfven speed, and the dotted line the sonic surface. 
Averaged over times t = 130. ..140. 
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Fig. 16. The ratio of the poloidal magneto-centrifugal and 

pressure forces, |-Fp™]^Vl-^poil defined in Eq. (p^, is 
shown with shades of grey, with larger values corresponding to 
lighter shades. The maximum value is l-F'p"^'^'' l/li^p^^l ~ 411. 
Superimposed are the poloidal magnetic field lines. The dashed 
line shows the fast magnetosonic surface with respect to the 
poloidal field (cf. Fig. The parameters are as in the model 
of Fig. ^ Averaged over times t = 130. ..140. 



celerated by pressure gradient near the disc surface (where the 
Alfven surface is close to the disc surface). However, acceler- 
ation remains efficient out to at least r ~ 1 within the conical 
shell at i9 « 30°. This can be seen in the upper and middle 
panels of Fig. [l^ (note that the conical shell is thinner and at 
a smaller i? in the model with larger entropy contrast, and so it 
cannot be seen in this figure, cf. Fig. These facts strongly in- 
dicate that magneto-centrifugal acceleration dominates within 
the conical shell. 

Another indicator of magneto-centrifugal acceleration in 
the conical shell is the distribution of angular momentum (see 
Fig. |l5|), which is significantly larger in the outer parts of the 
conical shell than in the disc, which suggests that the magnetic 
field plays an important role in the flow dynamics. We show in 
Fig. |l6|the ratio of the 'magneto-centrifugal force' to pressure 

gradient, l-fp^Y^ l/l-^poil' where the subscript 'pol' denotes the 
poloidal components. Here, the 'magneto-centrifugal force' in- 
cludes all terms in the poloidal equation of motion, except for 
the pressure gradient (but we ignore the viscous term and the 
mass production term, the latter being restricted to the disc 
only). 



g(n^z^-V<P) +JxB, 



(42) 



and i^*-^-* = — Vp. The large values of the ratio in the con- 
ical shell confirm that magneto-centrifugal driving is domi- 



nant there. On the other hand, the pressure gradient is strong 
enough in the outer parts of the disc to shift the Alfven sur- 
face close to the disc surface, leading to pressure driving. This 
is also discussed in Casse & Ferreira (2000b) who point out 
that, although the criterion of Blandford & Payne (1982) is ful- 
filled, thermal effects can be strong enough to lead to pressure 
driving. According to Ferreira (1997), a decrease of the total 
poloidal current /poi — i-KTuB^/ jii^ along a field line is an- 
other indicator of magneto-centrifugal acceleration. We have 
compared the poloidal current 1^^°^^ at the Alfven point, or at 
the top of our box if the Alfven point is outside the box, with 
the poloidal current /p^""^^' at the disc surface, and find that out- 
side the conical shell this ratio is typically > 0.8, while along 
the field line that leaves the box at (nj, z) — (0.6, 1), we get 



(top) I ^(surf ) 
pol 



/ po. 



0.18, i.e. a strong reduction, which confirms 



that magneto-centrifugal acceleration occurs inside the coni- 
cal shell. We note, however, that the changing sign of and 
therefore of /poi makes this analysis inapplicable in places, and 
the distribution of angular momentum (Fig. [I^ gives a much 
clearer picture. 

As further evidence of a significant contribution from 
magneto-centrifugal driving in the conical shell, we show in 



Fig. 17 that the magnetic field is close to a force-free configu- 



ration in regions where angular momentum is enhanced, i.e. in 
the conical shell and in the corona surrounding the outer parts 
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Fig. 17. The ratio J ■ B/\J\ \B\ in the corona for the model 
with Tstar = 0.01 and (i — 0.1, shown in Fig. 0. For a force- 
free magnetic field, J — CB, this ratio is ±1. Averaged over 
times t = 130. ..140. 




-0.5 



0.1 



figuration. The radial variation in the sign of the current helicity 
J ■ B is due to a variation in the sign of the azimuthal magnetic 
field and of the current density. Such changes originate from the 
disc where they imprint a corresponding variation in the sign of 
the angular momentum constant, see Eq. (p6|). These variations 
are then carried along magnetic lines into the corona. The loca- 
tions where the azimuthal field reverses are still relatively close 
to the axis, and there the azimuthal field relative to the poloidal 
field is weak compared to regions further away from the axis. 

Pressure driving is more important if the entropy contrast 
between the disc and the corona is larger (i.e. (3 is smaller): the 
white conical shell in Fig. indicative of stronger magneto- 
centrifugal driving shifts to larger heights for /3 — 0.02, as 
shown in Fig. |l8|. We note, however, that there are periods when 
magneto-centrifugal driving is dominant even in this model 
with higher entropy contrast, but pressure driving dominates 
in the time averaged picture (at least within our computational 
domain). 

We show in Fig. [l^ the variation of several quantities along 
a magnetic field line that has its footpoint at the disc surface 
at {zu,z) = (0.17,0.15) and lies around the conical shell. 
Since this is where magneto-centrifugal driving is still domi- 
nant, it is useful to compare Fig. |l^ with Fig. 3 of Ouyed et 
al. (1997), where a well-collimated magneto-centrifugal jet is 
studied. Since our outflow is collimated only weakly within 
our computational domain, the quantities are plotted against 
height z, rather than z/tij as in Ouyed et al. (1997) (z/w is 
nearly constant along a field line for weakly collimated flows, 
whereas approximately z/zu (x z along a magnetic line for 
well-collimated flows). The results are qualitatively similar, 
with the main difference that the fast magnetosonic surface in 
our model almost coincides with the Alfven surface in the re- 
gion around the conical shell where the outflow is highly super- 
sonic. Since we include finite diffusivity, the curves in Fig. [T^ 
are smoother than in Ouyed et al. (1997), who consider ideal 
MHD. A peculiar feature of the conical shell is that the flow at 
z ^5 1 is sub-Alfvenic but strongly supersonic. The fast mag- 
netosonic surface is where the poloidal velocity t;poi equals the 
fast magnetosonic speed for the direction parallel to the field 
lines. 



(43) 



1.0 



tO.Q 



IQO.Q 



with f A.poi the Alfven speed from the poloidal magnetic field. 
This surface has the same overall shape as the fast magne- 
tosonic surface with respect to the poloidal field, albeit in some 
cases it has a somewhat larger opening angle around the con- 
ical shell and is located further away from the disc in regions 
where the toroidal Alfven speed is enhanced. 



Fig. 18. As in Fig. |l6[ but for the model of Fig. ^ i.e. with 
larger entropy contrast, /3 = 0.02. Averaged over times t = 
140. ..240. 

of the disc. These are the regions where the Lorentz force con- 
tributes significantly to the flow dynamics, so that the magnetic 
field performs work and therefore relaxes to a force-free con- 



3.6. Lagrangian invariants 

Axisymmetric ideal magnetized outflows are governed by five 
Lagrangian invariants, the flux ratio, 

k = QU^JB^ = gu^/B^, (44) 
the angular velocity of magnetic field lines, 

h = w-\u^~kB^/g), (45) 
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Fig. 19. Flow parameters as functions of height z, along the 
field line with its footpoint at (-cufpjZfp) — (0.17,0.15), for 
the model of Fig. ^ Upper left: poloidal Alfven Mach num- 
ber. Ma = Wpoi/wA,poi (dashed), and a similar quantity that 
includes the poloidal Alfven speed as well as the sound speed 
(fast magnetosonic Mach number with respect to the poloidal 

)^/^ (dotted) - the two curves 



''A,pol 



field), Mfm = Mpoi/ (Cg 
are practically identical; upper right: ratio of toroidal (Stor) 
and poloidal (Spoi) magnetic fields; lower left: toroidal veloc- 
ity Mtor (solid) and poloidal velocity Upoi (dash-3dots), in units 
of the toroidal velocity at the footpoint, Utor.tp; lower right: 
density q in units of the density at the footpoint, gtp. The po- 
sition of the Alfven (and fast magnetosonic) point on the field 
line is indicated by the vertical line. This figure is useful to 
compare with Fig. 3 of Ouyed et al. (1997). 



the angular momentum constant, 

the Bernoulli constant, 

ft, + $ - wilB^/{^nk). 



(46) 



(47) 



and specific entropy s (which is a prescribed function of posi- 
tion in our model). In the steady state, these five quantities are 
conserved along field lines, but vary from one magnetic field 
Une to another (e.g., Pelletier & Pudritz 1992; Mestel 1999), 
i.e. depend on the magnetic flux within a magnetic flux sur- 
face, 27ra, where a = ruAip is the flux function whose contours 
represent poloidal field lines. 

We show in Fig. |20| scatter plots of k{a), f2(a), £{a), and 
e(a) for the model of Fig. ^ Points from the region 0.2 < z < 
8 collapse into a single line, confirming that the flow in the 
corona is nearly ideal.^ This is not surprising since the mag- 
netic Reynolds number is much greater than unity in the corona 

^ If we restrict ourselves to 0.2 < z 8, about 90% of the points 
for Q{a) are within ±10% of the line representing z = 4; for the other 
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Fig. 20. The four Lagrangian invariants fc(a), fl{a), £{a), and 
e(a) as a function of the flux function a for the model of Fig. ^ 
All points in the domain < ru < 8, 0.2 < z < 30 are 
shown (provided < a < 0.25). The dots deviating from well- 
defined curves mostly originate in 8 ^ z < 30. Note that spe- 
cific entropy as fifth Lagrangian invariant is trivially conserved 
in our model, because entropy is spatially constant throughout 
the corona and temporally fixed. 



for the parameters adopted here. For 8 SS z < 30, there are de- 
partures from perfect MHD; in particular, the angular velocity 
of magnetic field lines, ft, is somewhat decreased in the up- 



three invariants, this percentage is at least 97%. In a larger domain, 
0.2 < z ^ 30, these percentages drop to about 80% for k{a) and 
60% for Q,{a); for l{a) and e(a), however, they remain greater than 
90%. 
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per parts of the domain (indicated by the vertical scatter in the 
data points). This is plausibly due to the finite magnetic diffu- 
sivity which still allows matter to slightly lag behind the mag- 
netic field. As this lag accumulates along a stream line, the de- 
partures increase with height z. Since this is a 'secular' effect 
only, and accumulates with height, we locally still have little 
variation of fc and 17, which explains why magneto-centrifugal 
acceleration can operate quite efficiently. 

The corona in our model has (turbulent) magnetic diffu- 
sivity comparable to that in the disc. This is consistent with, 
e.g., Ouyed & Pudritz (1999) who argue that turbulence should 
be significant in coronae of accretion discs. Nevertheless, it 
turns out that ideal MHD is a reasonable approximation for the 
corona (see Fig. |o|), but not for the disc. Therefore, magnetic 
diffusivity is physically significant in the disc and insignificant 
in the corona (due to different velocity and length scales in- 
volved), as in most models of disc outflows (see Ferreira & 
Pelletier 1995 for a discussion). Thus, our model confirms this 
widely adopted idealization. 

4. Toward more realistic models 

The models presented so far have some clear deficiencies when 
compared with characteristic features relevant to protostellar 
discs. The first deficiency concerns the relative amount of mat- 
ter accreted onto the central object compared with what goes 
into the wind. A typical figure from the models presented above 
was that as much as 30% of all the matter released in the disc 
goes into the wind and only about 70% is accreted by the cen- 
tral object. Earlier estimates (e.g., Pelletier & Pudritz 1992) 
indicate that only about 10% of the matter joins the wind. 
Another possible deficiency is the fact that in the models pre- 
sented so far the low-temperature region extends all the way to 
the stellar surface whilst in reality the cool disc breaks up close 
to the star because of the stellar magnetic field. Finally, the 
overall temperature of the disc is generally too high compared 
with real protostellar discs which are known to have typical 
temperatures of about a few thousand Kelvin. 

The aim of this section is to assess the significance of vari- 
ous improvements in the model related to the above mentioned 
characteristics. We consider the effect of each of them sepa- 
rately by improving the model step by step. 

4.1. A larger sink 



As discussed in Sects. 3.1 and 3.2, properties of the outflow, 
especially of the nonmagnetic ones, are sensitive to the param- 
eters of the central sink. It is clear that a very efficient sink 
would completely inhibit the outflow near the axis. On the other 
hand, a magnetosphere of the central object can affect the sink 
efficiency by channelling the flow along the stellar magnetic 
field (Shu et al. 1994; see also Fendt & Elstner 1999, 2000). 
Simulating a magnetosphere turned out to be a difficult task 
and some preliminary attempts proved unsatisfactory. 

Instead, we have considered a model with a geometrically 
larger sink, and illustrate it in Fig. Here, rp — 0.25 in 
Eqs (|]) and (^3|), which is five times larger than the sink used in 
our reference model of Fig. ^. The relaxation time of the sink 
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Fig. 21. As in Fig. 0, but with a mass sink at the centre which 
is five times larger, ro = 0.25, and Tstar = 0.112. Note the 
absence of outflows along the axis and a larger opening angle of 
the conical shell, which crosses z — ±1 at tu w 1.2. Averaged 
over times t = 500 . . . 530, /3 = 0.1. 



in Eq. (jg), Tstar, was rescaled in proportion to the free-fall time 
at ro as Tstar ^ ro/vff oc rj]^^, which yields Tstar = 0.112. 
The resulting mass loss rate into the wind is Mwind ~ 0.8, 
which corresponds to 1.6 x 10~^Moyr~'^ in dimensional 
units. Although this is indeed smaller than the value for the ref- 
erence model (6 x 10~^ M© yr~^), the overall mass released 
from the disc is also smaller, resulting in a larger fraction of 
about 40% of matter that goes into the wind; only about 60% 
is accreted by the central object. As we show below, however, 
larger accretion fractions can more easily be achieved by mak- 
ing the disc cooler. We conclude that even a sink as large as 
almost twice the disc half-thickness does not destroy the out- 
flow outside the inner cone. However, the outflow along the 
axis is nearly completely suppressed. 

4.2. Introducing a gap between the disc and the 
sink 

In real accretion discs around protostars the disc terminates at 
some distance away from the star. It would therefore be unre- 
alistic to let the disc extend all the way to the centre. Dynamo 
action in all our models is restricted to vu > 0.2, and now we 
introduce an inner disc boundary for the region of lowered en- 
tropy as well. In Fig. ^ we present such a model where the 
Cdisc(f") profile is terminated at zu = 0.25. This inner disc ra- 
dius affects then not only the region of lowered entropy, but 
also the distributions of a, rjt, and g^"^'^. The radius of the 
mass sink was chosen to be ro = 0.15, i.e. equal to the disc 
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Fig. 22. As in Fig. ^ but with a gap between the sink and the 
disc with the disc inner radius at w = 0.25. The radius of the 
stellar mass sink is three times larger, 7-o — 0.15, and Tstar — 
0.052. Outflow is absent along field lines passing through the 
gap. Averaged over times t = 400 . . . 660, /3 = 0.1 in both the 
disc and the sink. 



semi-thickness. The correspondingly adjusted value of Tstar is 
0.052. The entropy in the sink is kept as low as that in the disc. 

As can be seen from Fig. |2[ the gap in the disc translates 
directly into a corresponding gap in the resulting outflow pat- 
tern. At the same time, however, only about 20% of the released 
disc material is accreted by the sink and the rest is ejected into 
the wind. A small fraction of the mass that accretes toward the 
sink reaches the region near the axis at some distance away 
from the origin and can then still be ejected as in the models 
without the gap. Note that our figures (e.g.. Fig. show the 
velocity rather than the azimuthally integrated mass flux den- 
sity (cf. Fig. |l3]); the relative magnitude of the latter is much 
smaller, and so the mass flux from the sink region is not signif- 
icant. 

4.3. A cooler disc 

We now turn to the discussion of the case with a cooler disc. 
Numerical constraints prevent us from making the entropy gra- 
dient between disc and corona too steep. Nevertheless, we were 
able to reduce the value of (3 down to 0.005, which is 20 times 
smaller than the value used for the reference model in Fig. 0. 
The (3 value for the sink was reduced to 0.02; smaller values 
proved numerically difficult. 

The value of /3 = 0.005 results in a disc temperature of 
3 X 10^ K in the outer parts and 1.5 x 10^ K in the inner 



parts. As in Sect. 4.2, the disc terminates at an inner radius 



Fig. 23. As in Fig. but with [3 — 0.005 in the disc and 
(3 = 0.02 in the sink. Averaged over times t — 150 . . . 280 
when the overall magnetic activity in the disc is relatively low. 



of vj = 0.25. At 137 « 0.5, the density in dimensional units is 
about 10~^ g cm""^, which is also the order of magnitude found 
for protostellar discs. The resulting magnetic field and outflow 
geometries are shown in Figs. 23 and ^ 

A characteristic feature of models with a cooler disc is a 
more vigorous temporal behaviour with prolonged episodes 
of reduced overall magnetic activity in the disc during which 
the Alfvenic surface is closer to the disc surface, and phases 
of enhanced magnetic activity where the Alfvenic surface has 
moved further away. Figures ^ and |4| are representative of 
these two states. It is notable that the structured outflow of the 
type seen in the reference model occurs in states with strong 
magnetic field and disappears during periods with weak mag- 
netic field. 

Another interesting property of the cooler discs is that now 
a smaller fraction of matter goes into the wind (10-20%), and 
80-90% is accreted by the central object, in a better agreement 
with the estimates of Pelletier & Pudritz (1992). 

We note in passing that channel flow solutions typical of 
two-dimensional simulations of the magneto-rotational insta- 
bility (Hawley & Balbus 1991) are generally absent in the sim- 
ulations presented here. This is because the magnetic field sat- 
urates at a level close to equipartition between magnetic and 
thermal energies. The vertical wavelength of the instability can 
then exceed the half-thickness of the disc. In some of our simu- 
lations, indications of channel flow behaviour still can be seen. 
An example is Fig. ^ where the magnetic energy is weak 
enough so that the magneto-rotational instability is not yet sup- 
pressed. 




Fig. 24. As in Fig. ^ but averaged over a period of enhanced 
magnetic activity in the disc, t — 300 . . . 320. A conical shell 
develops which crosses z — ±1 at tu sa 0.6. 

According to the Shakura-Sunyaev prescription, turbulent 
viscosity and magnetic diffusivity are reduced in a cooler disc 
because of the smaller sound speed (cf. Eq. (p^). Since in the 
corona the dominant contributions to the artificial advection 
viscosity fadv come from the poloidal velocity and poloidal 
Alfven speed (and not from the sound speed), c^'^'^ has to be 
reduced. Here we choose ass — 0.004 and reduce c^^^ by a 
factor of 10 to cff"^ = 0.002. Since we do not explicitly pa- 
rameterize the turbulent magnetic diffusivity % with the sound 
speed (cf. Eq. (^5])), also rjto needs to be decreased, together 
with the background diffusivity 770. We choose here values 
that are 25 times smaller compared to the previous runs, i.e. 
rjto = 4 X 10~^ and 770 = 2 x 10^^, which corresponds to agg^ 
ranging between 0.001 and 0.007. With this choice of coeffi- 
cients, the total viscosity and magnetic diffusivity in the disc 
have comparable values of a few times 10^''. 

The effect of reduced viscosity and magnetic diffusivity is 
shown in Fig. ^ Features characteristic of the channel flow 
solution are now present, because the vertical wavelength of 
the magneto-rotational instability is here less than the half- 
thickness of the disc. 

5. Discussion 

If the disc dynamo is sufficiently strong, our model develops a 
clearly structured outflow which is fast, cool and rarefied within 
a conical shell near the rotation axis where most of the angu- 
lar momentum and magnetic energy is carried, and is slower, 
hotter and denser in the region around the axis as well as in the 
outer parts of the domain. The slower outflow is driven mostly 



Fig. 25. As in Figs. || and ||, but with r/o ^ 2 x 10"^ 
T^to = 4 X 10-^ cf" = 0.002, and aa ^ -0.15. Averaged 
over times t = 230 . . . 236 (where the magnetic energy is en- 
hanced), ass = 0.004. 



by the entropy contrast between the disc and the corona, but the 
faster wind within the conical shell is mostly driven magneto- 
centrifugally. Without a central mass sink, the flow near the 
axis is faster, but otherwise the flow structure is similar to that 
with the sink. 

The half-opening angle of the cone with hot, dense gas 
around the axis is about 20°-30°; this quantity somewhat 
changes with model parameters but remains close to that range. 
The outflow in our models does not show any signs of col- 
limation. It should be noted, however, that not all outflows 
from protostellar discs are actually collimated, especially not 
at such small distance from the source. An example is the 
Becklin-Neugebauer/Kleinmann-Low (BN/KL) region in the 
Orion Nebula (Greenhill et al. 1998), which has a conical out- 
flow with a half-opening angle of 30° out to a distance of 25- 
60 AU from its origin. Therefore, collimation within a few AU 
(the size of our computational domain) is expected to be only 
weak. 

The region around the fast, cool and rarefied conical shell 
seen in Fig. ^ is similar to the flow structure reported by 
Krasnopolsky et al. (1999); see their Fig. 1. In their model, 
however, the thin axial jet was caused by an explicit injection 
of matter from the inner parts of the disc which was treated 
as a boundary. In our reference model the fast outflow is sub- 
Alfvenic because of the presence of a relatively strong poloidal 
field, whereas in Krasnopolsky et al. (1999) the outflow be- 
comes super-Alfvenic at smaller heights. Outside the conical 
shell the outflow is mainly pressure driven, even though the 
criterion of Blandford & Payne (1982) is fulfilled. However, as 
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Casse & Ferreira (2000b) pointed out, pressure driven outflows 
might dominate over centrifugally driven outflows if thermal 
effects are strong enough. 

In our model, matter is replenished in the resolved disc in 
a self -regulatory manner where and when needed. We believe 
that this is an improvement in comparison to the models of 
Ouyed & Pudritz (1997a, 1997b, 1999) and Ustyugova et al. 
(1995), where mass inflow is prescribed as a boundary condi- 
tion at the base of the corona. If we put q'^^^'^ = in Eqs ([l|) 
and (0), the disc mass soon drops to low values and the out- 
flow ceases. This is qualitatively the same behaviour as in the 
models of, e.g., Kudoh et al. (1998). 

We should stress the importance of finite magnetic diffusiv- 
ity in the disc: although poloidal velocity and poloidal magnetic 
field are well aligned in most of the corona, dynamo action in 
the disc is only possible in the presence of finite magnetic dif- 
fusivity, and the flow can enter the corona only by crossing 
magnetic field fines in the disc. 

An outflow occurs in the presence of both dipolar and 
quadrupolar type magnetic fields, even though fields with dipo- 
lar symmetry seem to be more efficient in magneto-centrifugal 
driving (cf. von Rekowski et al. 2000). The effects of the mag- 
netic parity on the outflow structure deserves further analysis. 

The dynamo active accretion disc drives a significant out- 
ward Poynting flux in our model. Assuming that this applies 
equally to protostellar and AGN discs, this result could be im- 
portant for understanding the origin of seed magnetic fields in 
galaxies and galaxy clusters; see Jafelice & Opher (1992) and 
Brandenburg (2000) for a discussion. We note, however, that 
the pressure of the intracluster gas may prevent the magnetized 
plasma from active galactic nuclei to spread over a significant 
volume (Goldshmidt & Rephaeli 1994). 

Our model can be improved in several respects. In many 
systems, both dynamo-generated and external magnetic fields 
may be present, so a more comprehensive model should in- 
clude both. We used an a^fl dynamo to parameterize mag- 
netic field generation in the disc because we have restricted 
ourselves to axisymmetric models. As argued by Brandenburg 
(1998), dynamo action of turbulence which is driven by the 
magneto-rotational instability can be roughly described as an 
a^^l dynamo. But this parameterization can be relaxed in three- 
dimensional simulations where one may expect that turbulence 
will be generated to drive dynamo action. Such simulations will 
be discussed elsewhere (von Rekowski et al. 2002). 

Since our model includes angular momentum transport by 
both viscous and magnetic stresses, it is natural that the ac- 
creted matter is eventually diverted into an outflow near the 
axis; this is further facilitated by our prescribed entropy gradi- 
ent at the disc surface. We believe that this picture is physically 
well motivated (Bell & Lucek 1995), with the only reserva- 
tion that we do not incorporate the (more complicated) physics 
of coronal heating and disc cooling, but rather parameterize it 
with a fixed entropy contrast. We include a mass sink at the 
centre which could have prevented the outflow, and indeed the 
sink strongly affects nonmagnetized outflows. We have shown, 
however, that the magnetic field can efficiently shield the sink 
and thereby support a vigorous disc wind. 



The assumption of a prescribed entropy distribution is a 
useful tool to control the size of the disc and to parameterize 
the heating of the disc corona. However, it should be relaxed 
as soon as the disc physics can be described more fully. The 
energy equation, possibly with radiation transfer, should be in- 
cluded. This would lead to a self-consistent entropy distribution 
and would admit the deposition of viscous and Ohmic heat in 
the outflow. In the simulations by von Rekowski et al. (2002), 
entropy is evolved. 

We believe that a mass source is a necessary feature of any 
model of this kind if one wishes to obtain a steady state. In 
the present paper the mass source is distributed throughout the 
whole disc to represent replenishment of matter from the mid- 
plane of the disc. Alternatively, a mass source could be located 
near to or on the domain boundary. 
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